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Abstract 

Dimensional regularization is incompatible with the standard covariant pro- 
jection methods that are used to calculate the short-distance coefficients in 
inclusive heavy quarkonium production and annihilation rates. A new method 
is developed that allows dimensional regularization to be used consistently to 
regularize the infrared and ultraviolet divergences that arise in these per- 
turbative calculations. We illustrate the method by calculating the leading 
color-octet terms and the leading color-singlet terms in the gluon fragmen- 
tation functions for arbitrary quarkonium states. We resolve a discrepancy 
between two previous calculations of the gluon fragmentation functions for 
the spin-triplet P-wave quarkonium states. 



1 



Typeset using REVT^ 



I. INTRODUCTION 



A rigorous theoretical framework for calculating inclusive production cross sections (and 
inclusive annihilation decay rates) for heavy quarkonium in QCD has recently been developed 
The cross section (or decay rate) is expressed as a sum of the product of short-distance 
coefficients and long-distance matrix elements. The short-distance coefficients can be calcu- 
lated by perturbation theory to any order in as{mc), the running coupling constant at the 
energy scale of the heavy quark mass. The long-distance factors are expressed as matrix 
elements in nonrelativistic QCD (NRQCD), an effective field theory that reproduces QCD 
to any desired order in the velocity of the heavy quark. The matrix elements scale in a 
definite way with v, the typical relative velocity of the heavy quark in quarkonium. The 
NRQCD factorization formalism therefore allows the inclusive production cross section to 
be calculated systematically as a double expansion in ^^(mc) and v. 

In calculating the short- distance coefficients beyond leading order in a^, ultraviolet diver- 
gences and infrared divergences inevitably arise and need to be regulated. In perturbative 
calculations, the most convenient method for regulating both ultraviolet and infrared diver- 
gences is dimensional regularization. The main advantage of this method is that it preserves 
many of the symmetries of a field theory, including gauge invariance and Lorentz invariance. 
It also often leads to simpler expressions for the finite parts of a perturbative calculation. 
One feature of dimensional regularization that is particularly convenient in a nonrenor- 
malizable effective field theory like NRQCD is that it automatically sets power ultraviolet 
divergences to zero. This makes it unnecessary to explicitly subtract these divergences as 
part of the renormalization procedure. Dimensional regularization has been widely used in 
perturbative calculations of the cross sections and decay rates of elementary particles. It has 
also been used in some perturbative calculations involving quarkonium, but the consistency 
of the method has not been carefully examined. 

Most calculations of production cross sections and annihilation decay rates for heavy 
quarkonium have been carried out using the covariant projection method [0. This method 
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involves the projection of a cc pair onto states with definite total angular momentum J. For 
orbital angular momentum L = 1 and higher, the projections are specific to 3 dimensions. In 
dimensional regularization, momentum integrals are analytically continued to = 3 — 2e di- 
mensions. There is therefore a potential inconsistency in combining the covariant projection 
method with dimensional regularization. Moreover, the projections for spin-singlet states, 
such as ^5*0 and ^Pi, involve the 75 matrix, whose extension to dimensions is problematic. 

We have recently developed an alternative method for calculating production cross sec- 
tions and annihilation decay rates for heavy quarkonium which fully exploits the NRQCD 
factorization framework 0. We refer to it as the threshold expansion method. In this method, 
a quantity that is closely related to the cross section for the production of a cc pair with 
total momentum P is calculated using perturbation theory in full QCD and expanded in 
powers of the relative 3-momentum q of the cc pair. Matrix elements of certain 4-fermion 
operators in NRQCD are also calculated using perturbation theory and expanded around 
the threshold q = 0. The short- distance coefficients in the factorization formula are then 
determined by matching these perturbative expressions. Finally, the NRQCD matrix ele- 
ments for quarkonium states are simplified by using rotational symmetry, heavy-quark spin 
symmetry, and the vacuum saturation approximation. 

In this paper, we generalize this method to N spatial dimensions, so that dimensional 
regularization can be consistently used to regularize infrared and ultraviolet divergences. The 
perturbative calculation of the QCD side of the matching condition and its expansion around 
the threshold generalize easily to dimensions. The perturbative calculation and threshold 
expansion of the NRQCD matrix elements can also be carried out in dimensions, provided 
that the matrix elements are defined in a way that generalizes to A^ dimensions. By matching 
the cross section and the matrix elements, we obtain the short-distance coefficients. After 
renormalization of coupling constants in QCD and NRQCD, the short-distance coefficients 
may have poles in A^ — 3, which must be removed by renormalization of the 4-fermion 
operators in NRQCD. One must take care to avoid simplifying the matrix elements of these 
operators using identities that are specific to 3 dimensions until after these renormalizations 
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have been carried out. 

In Section II, we review the threshold expansion method of Ref. [0], generalizing it 
to spatial dimensions. In Section III, we give a definition of the gluon fragmentation 
function for quarkonium states that is particularly convenient for low order calculations. 
We then present a matching prescription for calculating the short- distance coefficients in the 
factorization formula for these fragmentation functions. In Sections IV and V, we illustrate 
our calculational method by calculating the color-octet terms in the gluon fragmentation 
function for a general quarkonium state to order as and the color-singlet terms to order 
al- The color-singlet calculation involves a pole in — 3 which must be removed by the 
renormalization of an NRQCD matrix element. We apply our general formula for the gluon 
fragmentation function to S-wave and P-wave states in Section VI. We resolve a discrepancy 
between two previous calculations of the gluon fragmentation function for ^Pj states. 



II. THRESHOLD EXPANSION METHOD 

We begin by reviewing the threshold expansion method of Ref. [0], generalizing it to 
dimensions. Cartesian vectors and Lorentz vectors are extended to and A^ + 1 compo- 
nents, respectively. The inclusive cross section for producing the heavy quarkonium H with 
momentum P can be written in a factorized form : 

1 P 



where Ep = \J Mfj + P^. The coefficients Cmn take into account the effects of short distances 
of order 1/mc or smaller, and therefore can be calculated as perturbation series in the QCD 
coupling constant ^^(mc). The matrix elements (O^n) expectation values in the NRQCD 
vacuum of local 4-fermion operators that have the structure 

= X^I^U Vh ^^JCnX, (2) 

where ip and x ^^re the field operators for the heavy quark and antiquark in NRQCD, and 
JCn and JC'l^ are products of a color matrix (1 or T"), a spin matrix, and a polynomial in the 



gauge covariant derivative D in dimensions. The spin matrix is either the unit matrix 
or a polynomial in the Pauli matrices a*. The Pauli matrices in N dimensions satisfy the 
anticommutation relations 



{ (Ti, cTj } = 26ij, i,j = 1,...,N. (3) 

In 3 dimensions, they also satisfy the commutation relations 

[ai, Cj] = 2ieijk (Jk, i,j,k = l,2,3. (4) 

Using both the relations @ and the spin matrices in /C„ and )C'l^ in (H) can be reduced 
to a linear combination of 1 and a'. However, the commutation relations (D can be used to 
simplify the NRQCD matrix elements only after all poles in — 3 have been removed from 
the short-distance coefficients. Until these poles are removed, we must allow for additional 
spin matrices in )Cn and JC'l^, such as [cij, aj] and {[at, aj], ak}. The projection operator Vh 
in (|1|) can be written 

Vh = J2 |^(P = 0) + S){HiP = 0) + ^|, (5) 

s 

where the sum is over soft hadron states S whose total energy is less than the ultravio- 
let cutoff A of NRQCD. This operator projects onto the subspace of states which in the 
asymptotic feature include the quarkonium state H at rest plus soft hadrons. The standard 
relativistic normalization of the states in is 

H(P')\h(P)) = 2Ep {2TTf5^{P -P'). (6) 

With this normalization of states, the projection operator Vh has energy dimension 1 — A^. 

The short distance coefficients Cmn in ([1|) can be determined by matching perturbative 
calculations of the corresponding process in which the quarkonium H is replaced by a cc pair. 
Let cc(P, q, ^, rj) represent a state that consists of a c and a c with total momentum P, spatial 
momenta ±q in the cc rest frame, and spin and color states specified by the spinors ^ and 
rj. Color and spin indices on these spinors are suppressed. The number of spin components 



for these Pauli spinors is that appropriate to Pauh spinors in spatial dimensions. Using 
the abbreviated notation cc(P) = cc(P, q, ^,77) and cc'(P) = cc{P,q',C,',ri'), the matching 
condition in the threshold expansion method of Ref. is 



X 

/t 



pQCD 



= E Cmn {X^}CU 'Pc^,c, ^^^nX) (7) 

where ki and k2 are the momenta of the incoming particles and kx is the sum of the momenta 
of all the outgoing particles except the c and c. The operator Vcc',cc in the matrix element 
in (|^) is defined by 

Pcc-',cc- = E |c(q',r)c(-q',r/') + 5)(c(q, Oc(-q, r?) + S\. (8) 

s 

The sum is over soft parton states whose total energy is less than the ultraviolet cutoff A of 
NRQCD. The standard relativistic normalization is 

c(q;,0c(q^,^')c(qi,0c(q2,^)) = 4i?,,i?,, (27r)2^5^(qi - q;)<5^(q2 - q'^) ^W^r^, (9) 



where Eg = ym^ + q^. The spinors are normalized so that = 1, and similarly for t], ^' 
and 7]'. In expressions like both the spin and color indices are contracted. With the 
normalization (P), the projection operator Vcc',cc has dimension 2 — 2N. To carry out the 
matching procedure, the left side of (|^) is calculated using perturbation theory in full QCD, 
and then expanded in powers of q and q'. The matrix elements on the right side of (J^) are 
calculated using perturbation theory in NRQCD, and then expanded in powers of q and q'. 
The short-distance coefficients Cmn are obtained by matching the terms in the expansions 
in q and q' order by order in a^. 

In the perturbative calculations of the matching condition (|^, infrared and ultraviolet 
divergences can appear on both sides of the equation. Since NRQCD is constructed to be 
equivalent to full QCD at low energies, the infrared divegences on both sides must match. 
They therefore cancel in the short-distance coefficients Cmn- Any ultraviolet divergences 
on the left side are eliminated by renormalization of the QCD coupling constant and the 



heavy quark mass. On the right side, some of the ultraviolet divergences are eliminated 
by renormalization of the gauge coupling constant and other parameters in the NRQCD 
lagrangian. The remaining ultraviolet divergences are removed by renormalization of the 
4-fermion operators of NRQCD. 

The matching calculations are particularly simple if dimensional regularization is used to 
regulate both infrared and ultraviolet divergences. With dimensional regularization, power 
infrared diverences and power ultraviolet divergences are automatically set equal to zero. 
The only divergences which remain are logarithmic divergences, which appear as poles in 
e = (3 — A^)/2. With dimensional regularization, the NRQCD side of the matching condition 
is especially simple, because radiative corrections to the matrix elements vanish identically. 
The reason for this is that one must expand the integrand of the radiative correction using 
a nonrelativistic expansion in the loop momentum before integrating. Since the integrand 
is also expanded in powers of q and q', there is no momentum scale in the dimensionally 
regularized integral and it therefore vanishes. The radiative corrections have both infrared 
and ultraviolet divergences, but the infrared poles in e cancel the ultraviolet poles. Thus 
the only contributions to the NRQCD side of the matching condition are the tree-level 
contributions of the matrix elements, including those matrix elements that arise from coun- 
terterms associated with operator renormalization. The poles in e from the coefficients of 
the counterterm matrix elements will match the infrared poles in e on the QCD side of the 
matching condition. The short-distance coefficients of the renormalized matrix elements in 
the factorization formula depend on the operator renormalization scheme for NRQCD. The 
renormalized matrix elements in the minimal subtraction scheme are defined by the condi- 
tion that the coefficients of the counterterms in the expression for the bare operators in terms 
of renormalized operators are pure poles in e. In the more conventional MS renormalization 
scheme, the coefficients of the counterterms are multiples of l/e-\- log(47r) — 7, where 7 is 
Euler's constant. One must be wary of simplifying the matrix elements using identities that 
are specific to 3 dimensions until after these renormalizations have been carried out. 
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III. FACTORIZATION OF FRAGMENTATION FUNCTIONS 



To illustrate the use of dimensional regularization in the threshold expansion method, 
we will calculate the terms of lowest order in in the gluon fragmentation function for an 
arbitrary quarkonium state. The fragmentation function Di^niz) gives the probability for 
a jet initiated by a high energy parton i to include the hadron H carrying a fraction z of the 
jet momentum. The factorization formula ([I|), when applied to the fragmentation function 
for a gluon into a heavy quarkonium state H, yields the factorized form 

D,^Hiz) = J2dmn{z){O^J, (10) 

mn 

where the short- distance coefficients dmn{z) can be calculated using perturbation theory. 

Fragmentation functions can be calculated directly from the gauge-invariant field- 
theoretic definitions [Q. This method is particularly advantageous for higher order cal- 
culations. We choose instead to follow as closely as possible the calculational strategy 
introduced in Ref. 0, in which the fragmentation function is obtained by calculating the 
decay rate of a virtual parton in the infinite momentum frame and in an appropriate axial 
gauge. We consider the fragmentation of a gluon into a quarkonium state H. We take the 
gluon momentum to be £ = (£o, 0, . . . , 0, £jv). The infinite momentum frame is defined by 
the limit Iq^In oo with = — of order m^. We choose an axial gauge with reference 
vector n = (1, 0, . . . , 0, — 1), so that the gluon propagator is 

We also introduce n = (1, 0, . . . , 0, +1). Let "^g"_^j|/(p)+x be the amplitude for a virtual 
gluon to decay into a quarkonium state H plus additional final state particles X whose total 
momentum is kx- The gluon fragmentation function can be defined by 



-(A 

8{N - 1) (£2 



The sum over additional final state particles X includes integration over their phase space. 
The last factor in (|T2|) is the numerator of the gluon propagator in (|1^) evaluated at the 
positive energy pole d-o = d-N + i^- It can also be written in the form 

- I - N~l 

-^.. + "^"^^^^^'^ = E^M^ (13) 

i=l 

where e^^\ i = 1, . . . , N — 1, are the transverse polarization vectors for a real gluon whose 



momentum is proportional to n. The factor of 1/(8(A^ — 1)) in (0) comes from averaging 
over the color and polarization states of the decaying gluon. In calculations beyond leading 
order, there are infrared poles in A^ — 3 that arise from the splitting of the gluon into collinear 
partons. They should be absorbed into the Altarelli-Parisi evolution of the fragmentation 
functions. This complication will not enter into the leading-order calculations presented in 
this paper. 

The short-distance coefficients dmn{z) defined by the factorization formula (0) can be 
obtained by a matching prescription analogous to (0) for the cross section: 



27r J {2tt)^2Ep ^ V d-n 

1 If , n^Uy + n^n^ 



l{N-l) (£2)2 



dfiu H - I {-^3*^05' {P)+X J -^g* 



g'^cc(P)+X 



pQCD 



mn 

The left side is calculated using perturbative QCD and expanded around the thresholds 
q = q' = 0. The matrix elements on the right side are calculated using perturbative 
NRQCD and also expanded around the thresholds. Matching the expansions in q and q' 
order by order in as, we can determine the short-distance coefficients dmn{z). 



IV. COLOR-OCTET TERMS AT ORDER as 

The leading color-octet terms in the gluon fragmentation function Dg^H{z) have short- 
distance coefficients of order as- These terms come from the decay of the virtual gluon 
through the process g* —>■ cc. For this simple case, there are no additional final state particles 
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X in the left side of the matching condition (0). The energy-momentum-conserving delta 
function forces P = £ and the QCD side of the matching condition collapses to 

- --) 8(iv-'i)(n^ + "^"^2"'"^ (•^^•-<-))* -4r-(-.l«c. ■ (15) 

At leading order in a^, the amplitude for c(p)c(p) is given by the Feynman diagram 

in Fig. 1: 

^J"_,,(P) = gsf^'u{p)YT%{p). (16) 

The coupling constant in (^) has been written gsfJ-"^, where e = (3 — N)/2 and fi is the 
scale of dimensional regularization, so that Qs remains dimensionless in N dimensions. The 
momenta of the c and c are p = |P + Lq and p = |P — Lq, where is the matrix that 
boosts a spacelike vector in the rest frame of the cc pair to the frame in which the pair has 
total momentum P. The components of the boost matrix are given in of Appendix A. 



The spinor factor in (|T^) is expressed in terms of nonrelativistic Pauli spinors in ([A9b| ) . The 
amplitude becomes 

^r-cc-(P) = ^^c9sf^' ^a'T'^v ■ (17) 



Inserting this into the QCD side of the matching condition (|T5|), it reduces to 



Using the explicit expression for the boost matrices in ( |A4|) , the contraction of LjL'^ with 
Uufif^ reduces to {n ■ L)j{n ■ L)i = —z^z"^, where z is the unit vector in the A^'th coordinate 
direction. Using also the identity ( |A5a| ) , the QCD side of the matching condition reduces to 

5(1 - ^) g(J_;f)^, {5'^ - v'^a^rr^a^T^v. (19) 



We now consider the NRQCD side of the matching condition (p^. The spinor factor in 
(P!P|) can be identified as the expansion to leading order in as and to linear order in q and 
q' of the following NRQCD matrix element: 
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pNRQCD 



Ami r/'V^T'^^'^VT'^r/. 



(20) 



The tree-level expression for the matrix element in (^) is represented diagramatically in 
Fig. 2. The dot represents the operators ip'^a^T"'x X^^'^T'^'^i which create and annihilate 
cc pairs from the vacuum. The 2 lines emerging from the right side of the diagram represent 
the c and c in the bra of the projection operator Vcc',cc defined in (^. The 2 lines entering 
the left side of the diagram represent the c and c in the ket of the projection operator. The 
radiative corrections to the matrix element include diagrams in which virtual gluons are 
exchanged between the c and c on the right or between the c and c on the left. They also 
include diagrams in which real gluons enter the diagram on the left side and emerge on the 
right side with the same momenta and in the same color and spin states. These contributions 
come from higher Fock states in the projection operator Vcc',cc defined in (H). The sum over 
soft states S in the definition of Vcc',cc includes integrals over the momenta of the real gluons 
and sums over their spin and color quantum numbers. Since {x^cr^T°"il''Pcc',cc'ip^cr^T°'x) 
is a vacuum matrix element, the only diagrams that are allowed are those for which no 
propagators are cut by a vertical line through the dot. 

The coefficient of the factor ( pO]) in (|19|) is the short-distance coefficient for the matrix 
element. The same short-distance coefficient will hold for operators defined using the pro- 
jection Vh{\), where A specifies the polarization of the quarkonium state H. The color-octet 
term in the gluon fragmentation functions at leading order in is therefore 

,2e 



32(A^ - l)mi 



{5'^ - z'z^) {x^a^TyVnix) ^V^T^x). 



(21) 



The indices i and j range from 1 to N . We will find that, in order to obtain ultraviolet 
finite results at order a^, the matrix element in ([2T| ) will require renormalization. After that 



renormalization has been carried out, we can take the limit — 3. Our final result for the 
order-a^ term in the gluon fragmentation function is then 



Dg^H{\){z) = S{l-z) 



64m4 



(22) 
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The indices i and j now range from 1 to 3. We have set the scale of the running couphng 
constant in the short-distance coefficient equal to rric, since that coefficient is only sensitive 
to momenta on the order of rric or larger. The superscript (/i) on the matrix element is a 
reminder that it is a renormalized matrix element that may depend on a renormalization 
scale /i. For a specific quarkonium state, it may be possible to simplify the matrix element in 
(^) by using rotational symmetry, heavy-quark spin symmetry, and the vacuum saturation 
approximation. For example, if we sum over polarizations, we can use rotational symmetry 
to set 

A A 

For unpolarized quarkonium states, the calculation of the fragmentation function can be 
simplified by using rotational symmetry at an earlier stage of the calculation. Under the 
rotation group in dimensions, q and q' transform as vectors and ^, t], and t]' transform 
as spinors. Denoting the action of an element R of the rotation group by -R.q, R.C,, etc., the 
average over the rotation group of a function of these vectors and spinors can be defined by 

/(q,e,^,q',e,^') = JdR f{R.q, R.^, R.r], R.q', R.^', R.r]'), (24) 

where dR is the invariant integration element on the group, normalized so that / dR = 1. 
We can average the factor (-^g?^cc'(P)+x)*'^g*'-+cc(P)+x '^^ right side of ([T^) over rotations 
if we also average the projection operator Vcc',cc on the left side. Alternatively, we can leave 
Vcc'fic unchanged and simply restrict the sum on the left side to matrix elements for which 
{x'' IC'l^ipip'^ JCnX) is a scalar. 

Applying this rotational average to the QCD side of the matching condition in (0), it 
becomes 

^(1 - ^) si7-lW. ^ ^^^^^4^) ^ (25) 

The average over rotations of the spinor factor is 

'^p^JT^^^^ifHi = ^6'^ ri'^a^T'^i'^'^a^T^T]. (26) 
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Using the identity ( |A5b|) , the QCD side of the matching condition in (^) immediately 
collapses to 



2e 



5(1 - z) ri'^a^T^i'i'^a^T^ 



(27) 



The spinor factor is the lowest order expression for the matrix element 



pNRQCD 



Ami r^'V^T^C'^V^T^r?. 



(28) 



The coefficient of this term in ( PTl) is the short-distance coefficient for the matrix element. 
Thus the QCD side of the matching condition can be written 

6(1 -z) (xV^T> Vcc' cc ^^(y^T'^x) 



32iVm4 ' ' ^ "'-'"''^ 

The fragmentation function for an unpolarized quarkonium state H is therefore 

where the short-distance coefficient is 



(29) 



(30) 



5{l-z). 



(31) 



In the short-distance coefficient, we have set = 3 and set the scale of the running coupling 
constant equal to nic- The superscript (/i) on the matrix element in (|30D indicates that it is 
a renormalized matrix element. The projection operator Vh in the matrix element in (|30D 
is the sum over helicities of the projection operators defined in (^: Vh = J^x'PHiX}- The 
result ( pOD can also be obtained by using the identity (|23| ) in the unpolarized fragmentation 



function ([2^). 



V. LEADING COLOR-SINGLET TERMS AT ORDER al 



The leading color-singlet terms in the gluon fragmentation function Dg^ii{z) have short- 
distance coefficients of order a^. These terms come from the decay of the virtual gluon 
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through the process g* — > ccg. The phase space integrals and the energy-momentum- 
conserving delta function in the matching condition ( |T^ ) reduce to 



{27T)^2EpJ (27r)^2|k| 
1 /•! dz 



{271^+^ 6^+\i - P - k) 

die 



(32) 



2 Jo z{l - z) J (27r)^-i 'V z 1 ^ 

If we only require fragmentation functions that are summed over the polarizations of the 
quarkonium states, the calculation can be simplified by averaging both sides of the matching 
condition over rotations of the vectors and spinors that specify the state of the cc pair in its 
rest frame. The integrand on the QCD side of the matching condition ([l^) then reduces to 
a function of P^, and the integration over can be carried out using the remaining delta 
function in (|32D . The QCD side of the matching condition in (|14|) reduces to 



(v^)-(^+i) ds ^^_3 / n^n^ + n^n^ 



.2 ^' 



-9tiu + 



■g*^cc'{P)+g) •^g'^cc(P)+g 



pQCD 



, (33) 



16 r (^) ■/4m?/. S2 - - V ■ 2 

where s = ^ and P]_ = (1 — z){zs — 4Pg). At leading order in as, the amplitude for the 
process g*{l) c{p) c{p) g{k), with the cc pair in a color-singlet state, is the sum of the 
Feynman diagrams in Fig. 3: 
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e'^ik) u{p) 



vip). (34) 



2p ■ k 2p ■ k 

Using the identities for carrying out the nonrelativistic expansion of spinor factors given in 
Appendix A, we expand (|34D to linear order in q: 



'^g*^cc{P)+g 



6P-k 



+ 



rrir 



3 

Ami 



(k-L), [k>^Lr + L>'^r-{k-L)jg>^'' 



P-k 

- (k ■ LUP^'ir - L'^^P") + {k ■ ^.^P^^Lr - L'^.P'^) 

\^q'a^n]. (35) 



+ P-k L'^iL'^j 



P ■ i L^jL", 



In > 3 dimensions, the spin matrix {[a*, o"-'], cr^}, which is totally antisymmetric in its 
three indices, is linearly independent of 1 and a^,i = 1, . . . , N. In 3 dimensions, it reduces 
to the unit matrix multiplied by —Aie^^^. 
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After inserting the amplitude (|35D into the QCD side of the matching condition (p3D, 
the spinor factors can be simphfied by averaging over the rotation group. The factor 
^/tq,/fc^/^/^tg«(ji^ can be reduced to a hnear combination of three independent rotationally- 
invariant spinor factors: 

1 



^/f q'k^i ^/^t qi-(ji ^ 



N{N -1){N + 2) 



(36) 



The spinor factors ?7'''"g''cr™^'^''"{[(T*, a-'], cr'^}?7 and ?7''''{[cr', a™'], cr^j^'^^gV-^?] average to zero, 
while ?7'^{[cr*, cr-'], (t'^}^'^^{[ct*, cr-'], (j'^jr] reduces to a single rotationally-invariant spinor fac- 
tor: 



1 

N{N -1){N -2) 




5^' 






(37) 






^km 












^kn 







We will find that there are no poles in — 3 multipying this spinor factor. It can therefore 
be simplified by using the commutation relations (^) for A^ = 3: 

Vt{[a^a^],a*K'eUK,^1,^*h = -96VWr/, AT = 3. (38) 

After averaging the spinor factors over the rotation group using (^6]) and (^), all the 
Cartesian indices of the boost tensors Lf are contracted and they can be simplified using the 
identity ([A5b|) of Appendix A. The factor e'^'"(A;)*e'^"(/c), summed over colors and polariza- 
tions of the real gluon, can be replaced by 8{—g'^'''). After simplifying all the Lorentz algebra, 
the QCD side of the matching condition (^) reduces to integrals of scalar quantities: 



72 r (^) A^(A^ - 1)(A^ + 2) m4+2^ (1 + x)' 
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t 



(Af-3)/2 



where x = {s — 4m^)/(4m^), t = (1 — z)(za; + z — 1), and 



(39) 



W^o(x,^) 
W^i(a:,^) 



W2(X,Z) 



6{N + 2)x'^ -4t(l + a;) + (A^ - l)a;^ 



-At 



m + 4:{N - 2)x +{N - 9)x^ + {N - l)x'- 



+ {N- l)x^ 4 - 4(A^ + l)x + {N- l)x' 



(40a) 



(40b) 



-At 



-8 - 16x + {2N^ + N - 17)x^ + {2N^ + N - 9)x^ 



+ {N- l)x^ A{N^ - 3) + A{N^ - 3)x + {2N^ + N - 9)x^ 



(40c) 



W-,{x,z) = -At 



AN - A{N^ -2N - A)x - {6N^ - 5N - 23)x^ - {2N^ -N - 7)x'- 



+ {N- l)x^ A - A{N^ - 5)x - {2N^ -N- 7)x' 



(40d) 



Upon integrating ( |39D over x, we obtain distributions in z. Some of these distributions 
are infrared divergent in the sense that the integrals over z diverge as — > 3. The infrared 
divergences arise from the following integrals: 



00 f-c 

dx — 

il-z)/z 



r(l + e)r(l-e)^(l-^ 



dx 



t^-' r(2 + e)r(2-e) 



z' a 



\-l-2e 



i{i-z)/z X'* 6 ' ■ 

The divergences as e — > can be made explicit by using the expansion 



(41) 
(42) 



(1-^) 



-l-2e 



-6(1 -z) 
2e ^ ' 



1 



'l-z) 



(43) 



There is no infrared divergence in the Wq term. The divergences in the Wi and terms 
are multiplied by factors of e, and they give 5(1 — 2;) terms. There is an infrared divergence 
that survives in the W2 term. 

Integrating over x in ( P5| ) and taking the limit e — 0, the QCD side of the matching 
condition reduces to 
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+ d2{z,^) 



(44) 



The coefficient of the ffist spinor factor is 



at 



The remaining three functions in are 



[2(1 -z) \n{l - z) + z{3 - 2z)] . 



(45) 



di{z) 

C?2(^,/U) 



a: 



6480m6 

„2 



[5(1 - ^) + 6(17 - llz) ln(l - z) + 15^(7 - 2z)] , 



(46a) 



a: 



6480m6 



40 ■ 



1-z) 



40 In 



/i 



2m, 20 



36(2 - z) ln(l - 2) + 10z(5 - 4z) 



4 (-2) 



a: 



[5(1 - ;z) + 36(2 - z) ln(l - ;z) + 60z] . 



(46b) 
(46c) 



6480m6 

The subscript IR on the e in ( ^Ij ) is a reminder that the pole is of infrared origin. The 
reason for associated the term ln(47r/i^) — 7 + | with the pole in e in (^^) will become clear 
later. 

We now consider the NRQCD side of the matching condition (0). The spinor factors in 
(0) can be identified as the tree-level expressions for NRQCD matrix elements, expanded 
to linear order in q and q'. The ffist matrix element is 



pNRQCD 



Ami v'H'i'^V 



(47) 



and the remaining three are 



pNRQCD 
pNRQCD 
pNRQCD 
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4m^ V^q' ■ cr^'^^q ■ f^, 

4ml r^'^g'^a^e'^^g^fx"^?/ 



(48a) 
(48b) 
(48c) 



For the matrix elements (^T]), ( [48a|) , and ([48 c|) , we can immediately read off the short- 
distance coefficients from (0). They are the functions of z given in (P5|), ( [46a| ), and ( [46 c[ ). 
However the coefficient of v('^c('^o^^^^(f^o'^ri in (0) contains an infrared pole in e, indicating 
that it is sensitive to long-distance effects. Since an infrared divergence cannot appear in a 
short- distance coefficient, that divergence must be matched by an infrared divergence from 
some matrix element on the NRQCD side of the matching condition. Infrared divergences 
in NRQCD matrix elements can arise only from radiative corrections. Since the divergence 
in ( ^11) has a coefficient of order a^, the infrared- divergent NRQCD matrix element must 
have a short-distance coefficient of order as- The only such scalar matrix element is the 



one whose short-distance coefficient has already been determined in (p9|). Thus the infrared 
divergence on the NRQCD side of the matching condition must come from that term. 

If the cc pairs are in color-singlet states, the tree level expression ( p8[ ) for the matrix 
element vanishes and the leading contribution comes instead from radiative corrections. 
Specifically, it comes from the ccg term in the projection operator Vcc',cc defined in (||): 

E J ^^^|c(q^m-q^V)+^7(k,A,6))(c(q,Oc(-q,r/) + <7(k,A,6)|. (49) 

The sum is over the — 1 physical polarizations and the 8 color states of the real gluon. 
The leading contributions to the matrix element are represented by the diagrams in Figure 
4. The expression for diagram 4a is 

r h 

4//.^^ r^'^a-T'^n'^a-T'T'^r^ j (q ■ q' - q ■ k k ■ q') 

1 1 

(50) 



Eg + k-{q + k)V(2m,) + Eg + k - {q' + k)V(2m,) + ' 

where Eg = q^/(2mc) = (q')^/(2mc). As discussed in Appendix B of Ref. [|], the proper 
way to evaluate the diagram is to first expand out the denominators in powers of q/m^, 
q' /rric, and k/mc, and then integrate over k. Keeping only terms up to linear order in q/rric 
and q'/rric, the diagram reduces to 

C^^ , 2e ^nrriarrib c' C'i ^nrriarrib^ f ^ q " q' — q " k k ■ q' , . 

87ra,/i ri^a T T ^^^a T T 7] J j^-y^ — . (51) 



The integral is both ultraviolet and infrared divergent. It vanishes in dimensional regular- 
ization due to a cancellation between an ultraviolet pole in e and an infrared pole. Making 
these poles explicit, the diagram can be written 

-) r^'^q"^a''T^T'i'i'^q^a''T'T''7] . (52) 

3vr \euv (^irJ 

The subscripts UV and IR on e indicate whether the pole is of ultraviolet or infrared 
origin. We have set = 3 in the prefactor, since any finite terms obtained by expanding 
the prefactor in powers of e will cancel. The effect of the other 3 diagrams is simply to 
symmetrize both of the products of color matrices T°'T^. Since the projector Vcc',cc requires 
the asymptotic cc pairs to be in color-singlet states, we can replace {T",T^} by (5"^/3. The 
final result for the matrix element is 

(xV'TVP„,,,^V'T.x)t„„,, = W (^ - J^) n'mV'T",. (63) 

After multiplying by the short-distance coefficient in (^9]) , we find that the infrared pole in 
e matches the one on the QCD side of the matching condition, which is given in (^i]) . 

After taking into account the matrix element (^) on the NRQCD side of the matching 
condition, the short- distance coefficient of the matrix element ( [48b|) can be read off from 
(0). The net effect of taking into account the matrix element ( p3D is simply to change 
the infrared pole in e into an ultraviolet pole. The short-distance coefficient of the matrix 
element ( [48 b| ) therefore contains an ultraviolet divergence. This divergence must be removed 
by operator renormalization in NRQCD. As is evident from the ultraviolet pole in e in (^3|), 
it is the matrix element {x^ a^T^-ipVcc' ,cc4''^ <^^T°'x) that requires renormalization. In the MS 
renormalization scheme, the relation betweeen the matrix element of the bare operator and 
matrix elements of renormalized operators is 

+ + ln(47r) - 7) (x^(-fB)™a>Pec-',cEV'^(-tS)-a"x)^'^^) ■ (54) 

The superscripts (/i) on the matrix elements on the right side indicate that they are renor- 
malized matrix elements with renormalization scale /i. We will suppress this superscript on 
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P-wave matrix elements, since they do not require any renormalization at this order in a^. 
The fermion field operators in the bare matrix element on the left side of (^4|) have dimen- 
sion N/2. The fermion field operators in the renormalized matrix elements on the right side 
have dimension 3/2. The factor of on the right side of ( |5^ compensates for the differ- 
ence between the dimensions of the two sides. Solving (jS^) for {x^ a'^T'^ipVcc' ,cc'4''^ <^''T°'x)^^^ 
and using (|53|) and (|48b|) , we find that the renormalized matrix element, with dimensional 



regularization as the infrared cutoff, is 



(m) 



pNRQCD 



277T KeiR 



Multiplying ( pi]) by the short- distance coefficient in (29), we find that the contribution from 



the renormalized matrix element to the NRQCD side of the matching condition is 



32A^m4 



c 



pNRQCD 



"° ^ ^ +\n{Anfi'')-^ + l)6{l-z)Aml7]'^q'"'a''^'^^qV7]. (56) 



324m6 \euv 3. 

This term on the NRQCD side of the matching condition matches the part of the 
ri''^q'"^a"'C,'C,^q'^o'^Ti term on the QCD side that contains the pole in e. Subtracting this 
term from (^^, we find that the short-distance coefficient of the matrix element ( [48b|) is the 
function d2{z,^) given in ( [46b|) . 



Note that the calculation of the radiative correction to the NRQCD matrix element 
in (^) was necessary only to determine the coefficient of the counterterm in (|5^). The 
radiative correction itself vanishes if we identify euv = ^ir- Its only effect is to transform 
the ultraviolet pole in the coefficient of the counterterm into an infrared pole on the NRQCD 
side of the matching condition. 

Combining the color-octet term in (^) with the color-singlet terms determined above, we 
obtain a general expression for the gluon fragmentation function of an unpolarized quarko- 
nium state H: 
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where the short- distance coefficients are given in (pTD, (^51) , and (^q) . This expression 
for the fragmentation function includes all those color-octet matrix elements with short- 
distance coefficient of order and all those color-singlet matrix elements with short- distance 
coefficients of order a1 that involve at most one derivative acting on and ip and at most 
one derivative acting on and x- 

The /i-dependence of the coefficient d2{z,fi) in (|57| ) is cancelled by the /i-dependence of 
the color-octet matrix element. The renormalization group equation that determines the 
/i-dependence of that matrix element is obtained by differentiating ( p^ ) with respect to fi. 
We use the fact that the bare matrix element on the left side is independent of fi and that 
the bare coupling constant, which at this order is asfi'^'^, is also independent of fi. Replacing 
Vcc',cc by the projector Vh, we find 



^ij^ix^a'^TyVH^^cr'T'^x)^^^ = -^a,ifi){xK-0ra'^^PVHijK-0rcT'^x)- (58) 



The solution to this renormalization group equation at leading order in is 



{x^(r''T^iljVHip^(r''T''xY^^ = (x^^^'^> 



32 log f^l^\ (^^1^_0ym^n^ ^t(_iB)-^n^) , (59) 



9(33 — 2n/)m^ \as(/io 
where ri/ is the number of flavors of light quarks (n/ = 3 for charmonium). 

VI. GLUON FRAGMENTATION FUNCTIONS 

The relative importance of the various terms in the fragmentation function (|57D depends 
on the quarkonium state. The magnitude of a particular term is determined by the order in 
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as of its short- distance coefficient and by the scahng of the matrix element with v, which is 
given by the velocity-scahng rules of NRQCD [Q. Below, we apply this general fragmentation 
formula to specific S-wave states and P-wave states. 

A. Spin-singlet S-wave states 

The dominant Fock state of the rjc consists of a cc pair in a color-singlet ^5*0 state. 
The dominant matrix element is therefore (x^'ipVr^^ip'^x) ^ which scales as v^. It has a short- 
distance coefficient of order a^, and therefore contributes to the gluon fragmentation function 
at order a^v^. Since the dominant Fock state can be reached from a color-octet ^5*1 state 
through a chromomagnetic dipole transition, the matrix element {x^ a^T'^'ipVr^^'ip'^ a^T"'x) is 
suppressed by relative to the dominant matrix element. It contributes to Dg^^{z) at 
order UsV^. The color-singlet P-wave matrix elements in (|5^) are all suppressed by f ^. Their 
contributions are therefore of order a^f and are expected to be negligible. 

The most important term in the gluon fragmentation function for the r/c is the color- 
singlet ^5*0 term. Keeping only this term, the fragmentation function reduces to 

(60) 

Since the short- distance coefficient is dominated by contributions from large momenta of 
order or larger, we have set the scale of the running coupling constant equal to m^. 
The standard NRQCD matrix elements {Oni^^^^Lj)) introduced in Ref. [|l| were defined 
using a projection operator analogous to Vh in (H), except that the states have the standard 
nonrelativistic normalization. The relation between our matrix elements and those defined 
in Ref. is discussed in Appendix B of p. At leading order in f^, they differ simply by a 
normalization factor. For the matrix element in (|60D, the relation is 

(xV^e^^x) ^ Amnion's,)). (61) 

The term ( |5DD in the fragmentation function was first calculated by Braaten and Yuan [Q. 
The next most important term should be the color-octet ^Si term in (|57|). 
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D 



3z - 2^2 + 2(1 - z) ln(l 



B. Spin-triplet S-wave states 



The dominant Fock state of the J /ip consists of a cc pair in a color-singlet ^5*1 state. The 
dominant matrix element is therefore {x^a^ipV^ip'^a^x)) which scales as v^. Since the domi- 
nant Fock state can be reached from a color-octet ^5*1 state through a double chromoelectric 
dipole transition, the matrix element [x^ a^T""il)V^ip'^ a^T°'x) is suppressed by relative to 
the dominant matrix element. The matrix element (xV^^/'V'^x) is suppressed by f and the 
color-singlet P-wave matrix elements in (^) are all suppressed by v^. Their contributions 
are therefore expected to be negligible. 

Of the terms that we have calculated, the most important is the color-octet ^Si term. 
Keeping only this term, the fragmentation function reduces to 

D,_^^{z) = 5{1 - z)^^^^ {x^a^T^^V^^^a^T'^x). (62) 

Since the short-distance coefficient is dominated by contributions from large momenta of 
order iric or larger, we have set the scale of the running coupling constant equal to iric. 
Up to corrections of relative order f^, the matrix element in ( |62D differs from the standard 
NRQCD matrix element introduced in Ref. [|l| only by a normalization factor: 

{x^a'T''^V^^^a'T\) ^ 4m, {OtCS,)) . (63) 

The result ( |62D was first given by Braaten and Fleming 0. This contribution to the frag- 
mentation function is of order a^f^. If we assume that the effect of a suppression factor 
of as is comparable to that of a suppression factor of f^, then there is one other term in 
the fragmentation function that could be equally important. The color-singlet ^5*1 term 
contributes at order a^t>^. It has been calculated in Ref. 0. 

C. Spin-singlet P-wave states 

The dominant Fock state of the he consists of a cc pair in a color-singlet 
^Pi state. According to the velocity-scaling rules, the largest matrix elements are 
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(x"^(-|D)™V^fe^V^t(_|D)™x) and {x^T''i:Vh,i^^T''x) , both of which scale as v^. The rea- 
son that they are comparable in magnitude is that one is suppressed by a factor of f ^ from 
the covariant derivatives, while the other is suppressed by f ^ because of the chromoelectric 
dipole transition required to reach the dominant Fock state from a color-octet ^5*0 state. 
The color-octet ^Si and color-singlet ^5*0 matrix elements in (|57D contribute at orders asV^ 



and a^v^, respectively. The color-singlet P-wave matrix elements in (^) contribute at order 
and should be negligible. If we assume that the effect of a suppression factor of as 
is comparable to that of a suppression factor of f^, then the most important term is the 
color-octet ^5*0 term, which contributes at order a'^v^. This term has not yet been calculated. 

D. Spin-triplet P-wave states 

The dominant Fock state of the Xcj consists of a cc pair in a color-singlet ^Pj state. Ac- 
cording to the velocity-scaling rules, the largest matrix elements are {x^ a^T"-tpV^^jip^ cr^T"'x) 
and the color-singlet P-wave matrix elements in (|57|) , all of which scale as v^. The color- 
octet ^5*1 matrix element contributes to the gluon fragmentation function at order agV^, 
while the color-singlet P-wave matrix elements contribute at order a^v^. The matrix ele- 
ment {x^ipT^xcj'^'^X) contributes at order a^f ^ and should be negligible. We will therefore 
consider only the color-octet ^Si term and the color-singlet P-wave terms in the fragmenta- 
tion function. 

The symmetries of NRQCD can be used to reduce the three color-octet S-wave matrix 
elements {x^ a'^T^'ipV^^jip'^ a'^T'^x) to one independent matrix element, which we take to be 
the J = one. These symmetries are rotational symmetry and heavy-quark spin symmetry. 
Spin symmetry is an approximate symmetry with corrections that are of relative order v^. 
It implies that the state Xcj(A), which is an eigenstate of and J^, can be expressed in the 
form 

Xcj(A)) ^ ^(l/,;ls,|JA)|xc(/.s.)), (64) 

IzSz 
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where xMzSz) is an eigenstate of and Sz- Using the definition (|) of the projection 
operator Vh, the color-singlet ^5*1 matrix element can be written 

X {0\x^cT'T^^\xcihSz) + S){xc{l'A) + S\ij^a'T^x\^) . (65) 
s 

Spin symmetry also implies that the two matrix elements in ( |65D are proportional to Us^i 
and f//^, , where Umi is the unitary 3x3 matrix that transforms vectors from the Cartesian 
basis to the spherical basis. Finally, rotational symmetry implies that the product of the two 
matrix elements in ( |B5D summed over soft states S must be proportional to Sij>^. Using the 
orthogonality relations of the Clebsch-Gordan coefficients, the equation (|65|) can be reduced 
to 

ix^a'Ty V^^, ^ VT'^x) ^ (2 J + 1) (xV'T> V^^, ^^a'T^x) ■ (66) 

These relations hold up to corrections of relative order f ^. 

The color-singlet P-wave matrix elements in (|57D can also be reduced to a single inde- 
pendent matrix element, which we choose to be fD • (T)ipV^^f^ip'' {—^D ■ cr)x)- Using 

the expression (Q) for the Xcj states and the definition (R), we can write 



X E (0|x^(-f D)V^-^|xc(45.) + S){xc{l'-A) + D)-a"x|0) • (67) 

s 

Spin symmetry implies that the two matrix elements in (^) are proportional to Ug^j and 
Ul^gi . Since the operators create and annihilate cc pairs in the dominant Fock state of Xc, we 
can use the vacuum-saturation approximation, which is accurate up to corrections of relative 
order v'^. Keeping only the vacuum term in the sum over 5, the matrix elements reduce to 
(0|x"^(-iD)V^V'|Xc(4s^)) and (xc(/^4)|^^(-|D)™(t"x|0). By rotational symmetry, these 
must be proportional to Ui^i and f/^^/ . Thus the tensorial structure of the matrix element 
(|57|) is completely determined. The proportionality constant can be deduced by taking the 
special case J = 0, i = j and m = n. The resulting formula is 
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X E E E(l^- ls,\JX){JX\K; Is'M.jUlUi^^Ul, . (68) 

The scalar combinations of these matrix elements can be simplified by using the orthogo- 
nality relations of Clebsch-Gordan coefficients together with the identity 



(UU^) = -V3 (lmi;lm2|00). (69) 

The resulting formulas are 

(xn-fD ■ cT)^V^^,i^K-\Y) . ct)x) ^ 6jo {x\--fi ■ ct)^IjV^^J^ {-^B ■ ct)x) , (70a) 
(xn-tD)-a>P^^>t(_iD)™a"x) ^ (x^-fD ■ ^)^P,>t(_iD . ^)^) ^ (70b) 

(xn-tD)-cr>P^^>t(_iD)V-x) 

^ (-1)^^^ (x^-f D ■ a)^V^^,^\-^fi ■ a)x) . (70c) 

Using the relations ( pBD and ([7D|), the dominant terms in the gluon fragmentation function 
for the Xcj reduce to 

Dg^^^^,{z) = (2J+ 1) d^'-'^\z) {x^a-T^^ V^^, ^^a^T^x)^^^ 

+ ci(i''^-)(z,/x) (x^(-fD ■ ct)^ V^, ^t(_iD . ct)x) , (71) 

where S-'^^'^^z) is given in (|31|) and S^'^j\z,fi) is a linear combination of the functions in 



oTii 27-1-1 
d(i''^-)(z,/i) = 5jorfi(^) + -^d2iz,fi) + (-ly^^dsiz). (72) 



More explicitly, these coefficients are 

(2J+1) 



^ '^^ 486m6 



(l-^)-f 

+ (^Qj-(2J+l)log^)5(l-^) + Pj(z) 



(73) 



where the numbers Qj are 
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Qo = I, Qi = I, Q2 = I, (74) 

and the functions Pj{z) are 

^ , , ^(85 - 26z) 9(5 - 3z) , , 

Po{z) = ^ + ^ log(l - z) , (75) 

3^(1 + 42:) 

^"1(2:) = ^ , (76) 

P,{z) = + 9(2 - log(l - z) . (77) 

The matrix elements in (|7lD are related to the standard matrix elements defined in Ref. 
by 

(xVT> V^, V^VT^x) ^ 4m, Of (^Si) , (78) 
(x^(-f D ■ o-)^ ^t(_|D . ct)x) ^ 12m, Of (=^Po) • (79) 

We now compare our results for the gluon fragmentation functions of the Xcj states 
with previous calculations. These fragmentation functions were first calculated by Braaten 
and Yuan They regularized the infrared divergence from the process g* ccg by 
imposing a cut-off |k| > A on the momentum of the real gluon in the final state, and 
they used the covariant projection method to isolate the contributions to the fragmentation 
functions for Xco, Xa, and Xc2- They also assumed implicitly that the renormalization of the 
NRQCD matrix element {x^a'T''^Vcc',ccip^(r'T''x) 

was carried out in such a way that the 
net effect of terms on the NRQCD side of the matching condition was simply to transform 
the infrared cutoff A into an ultraviolet cutoff on NRQCD. The fragmentation functions 
calculated in Ref. [§] have been checked by Cho, Wise, and Trivedi p, who extended the 
calculation to the polarized fragmentation functions of g xj- Our results differ from those 
in Ref. p| only in the coefficients of the 6{1 — z) terms in ([731). Braaten and Yuan obtained 



Q'j - (2J + 1) ln(A/me), where Q'q = ^, Q\ = f , and Q'2 = ^. The numbers Qj and Q'j 
differ by |(2J + 1). Since the difference between the coefficients of 5(1 — z) is proportional 
to 2J + 1, it can be absorbed into a shift in the value of the color-octet matrix element 
{x^^ a^T^'ipV^gip'^ a'^T"'x) ■ The relation between the matrix elements in the two calculations is 
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dim.reg. 



cutoff 



Thus our calculation is consistent with that of Braaten and Yuan. The differences in the 
short- distance coefficients are due simply to different definitions of the color-octet matrix 
element. 

The gluon fragmentation functions for Xcj have also been calculated by Ma [|l^ . Ma used 
dimensional regularization to cut off the infrared divergence from the process g* — >■ ccg. He 
also used the covariant projection method to isolate the contributions to the fragmentation 
functions for Xco, Xci, and Xc2- His results differ from ours only in the numbers Qj multiply- 
ing the 5(1 — z) terms in (|73|). He obtained Qo = Qi = and Q2 = The differences 
between his values of Qj and ours are not proportional to 2J + 1. Therefore, they cannot be 
absorbed into a redefinition of the color-octet matrix element. We attribute the discrepancy 
to an inconsistency in Ma's calculation. Dimensional regularization of the infrared diver- 
gence involves analytically continuing integrals to 3 — 2e dimensions, while the standard 
covariant projection method uses projection matrices that are specific to 3 dimensions. 



VII. CONCLUSIONS 

Dimensional regularization is the most convenient method for regularizing the infrared 
and ultraviolet divergences that arise in calculations of quarkonium production and annihila- 
tion rates beyond leading order. The standard covariant projection method is incompatible 
with dimensional regularization. An alternative method for calculating the short- distance 
coefficients in the NRQCD factorization formulas is the threshold expansion method devel- 
oped in Ref. 0]. In this paper, we have generalized this method to N spatial dimensions, 
so that dimensional regularization can be used consistently to regularize the infrared and 
ultraviolet divergences that arise in perturbative calculations. We illustrated the method by 
calculating the color-octet terms of order as and the color-singlet terms of order in the 
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gluon fragmentation functions for arbitrary quarkonium states. We resolved a discrepancy 
between two previous calculations of the gluon fragmentation functions for the spin-triplet 
P-wave quarkonium states. This general and systematic method should be very useful in 
extending calculations of the production and decay rates of heavy quarkonium states to 
higher orders in perturbation theory. 
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APPENDIX A: NONRELATIVISTIC EXPANSION OF SPINORS 

This appendix is identical to Appendix A of Ref. except that it includes only those 
formulas that generalize to spacial dimensions. Formulas involving Levi-Civita tensors 
that cannot be easily generalized are omitted. We give the nonrelativistic expansions for the 
spinors of a heavy quark c and antiquark c for arbitrary momentum P = {Pq, Pi, ... , P/v). 
We assume that the relative momentum q = (gi, . . . , q^) of the c in the center-of-momentum 
(CM) frame of the cc pair is small compared to the quark mass rric. The momenta p and p 
of the c and c can be written 

p = |P + Lq, (Ala) 
p = iP - Lq , (Alb) 

where P is the total momentum and L is a Lorentz boost matrix. From the mass-shell 



conditions, p^ = p^ = m^, we have P ■ Lq = and P^ = ^E"^, where = \Jm1 + q^. The 
components of the momenta P and Lq in the CM frame of the pair are 



p/j. 
(Lq)^ 



= (2^„ 0) , (A2a) 

CM 

= (0, q) . (A2b) 

CM 
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When boosted to an arbitrary frame in which the pair has total spacial momentum P, these 
momenta are 



P>^ = [^^El + P2, P) , 



(A3a) 
(A3b) 



The boost matrix L'*-, which has one Lorentz index and one Cartesian index, has components 



2E. 



pi pj pO pi pj 
+ 

P2 ^ 2E„ P2 



(A4a) 
(A4b) 



The contraction of the boost tensor L'^j with the Lorentz vector P vanishes: P^L^j — 0. The 
contractions of two boost matrices in their Lorentz indices or in their Cartesian indices have 
simple forms: 



ptipv 



(A5a) 
(A5b) 



The representation for gamma matrices that is most convenient for carrying out the 
nonrelativistic expansion of a spinor is the Dirac representation: 



^ 0^ 



7 



-1 



^ A 



-a' 



(A6) 



In the CM frame of the cc pair, the spinors for the c and the c are 

1 



u{p) 



CM 



CM 



rrir 



\ 



rrir 



( \ 
— q ■ cr 77 

^ {Eq + mc)rj ^ 



(A7a) 



(A7b) 



Color and spin quantum numbers on the Dirac spinors and on the Pauli spinors ^ and rj are 
suppressed. When boosted to a frame in which the pair has total spacial momentum P, the 
spinors for the c and c are 
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u{p) 



1 



AE,{P^ + 2E^){E, + m, 
1 

AE,{P^ + 2E,){E, + m, 



{2Eg +f 7o) 



{2E^ +f 7o) 



V 



— q ■ (T T] 

^ {Eg + mc)ri ^ 



(A8a) 



(A8b) 



These spinors are normalized so that uu = —vv = 2mc if the Pauh spinors are normahzed 
so that = rj'^T] = 1. 

The independent quantities that can be formed by sandwiching 3 or fewer Dirac matrices 
between u{p) and v{p) are 



u 



u 



{p)v{p) = -2 ^■^(q- cr)r7, 
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The simplest way to derive these formulas is to use the identities 



7] . 
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Using the expressions for the spinor factors given in ([X9|), it is easy to carry out the nonrel- 
ativistic expansions in powers of q. 
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FIGURES 



Fig. 1. The lowest-order Feynman diagram in QCD for g*{l) — > c{p) + c{p). 

Fig. 2. The lowest-order Feynman diagram for the NRQCD matrix element 
{x^ a^T°-'4>Vcc' ,cc^'^ (^^T"''X) ■ The hnes entering at the left and leaving at the right come from 
a term of the form |cc')(cc| in the projection operator Vcc',cc- 

Fig. 3. The lowest-order Feynman diagrams in QCD for g*{l) — > c(p) -\-c{p) + g{k) when 
the cc pair is in a color-singlet state. 

Fig. 4. The lowest-order Feynman diagrams for the NRQCD matrix element 
{x^ a^T"-'4>Vcc' ,cc4''^ <^^T°'\) when the projection operator Vcc'^c requires the cc pairs to be 
in color-singlet states. The lines entering at the left and leaving at the right come from a 
term of the form |cc'g')(ccg'| in the projection operator. 
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